package it.unibo.refolding.pico;

import it.unibo.refolding.common.CommonThings;
import it.unibo.refolding.common.MyRuntimeException;
import it.unibo.refolding.pico.PicoDataAnalyser.ContourLengthEstimator;

import java.util.logging.Level;
import java.util.logging.Logger;

import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
import org.apache.commons.math3.analysis.MultivariateVectorFunction;
import org.apache.commons.math3.exception.TooManyEvaluationsException;
import org.apache.commons.math3.optim.InitialGuess;
import org.apache.commons.math3.optim.MaxEval;
import org.apache.commons.math3.optim.PointVectorValuePair;
import org.apache.commons.math3.optim.nonlinear.vector.ModelFunction;
import org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian;
import org.apache.commons.math3.optim.nonlinear.vector.Target;
import org.apache.commons.math3.optim.nonlinear.vector.Weight;
import org.apache.commons.math3.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer;

import com.google.common.base.Preconditions;

class LevenbergMarquardtOptimizerContourLengthEstimator implements ContourLengthEstimator {
  private final double kbtOverP;
  private final LevenbergMarquardtOptimizer estimator;

  private static final Logger logger = Logger.getLogger(LevenbergMarquardtOptimizerContourLengthEstimator.class.getSimpleName());

  public LevenbergMarquardtOptimizerContourLengthEstimator(double kbtOverP) {
    this.kbtOverP = kbtOverP;
    estimator = new LevenbergMarquardtOptimizer();
  }

  @Override
  public double estimate(final double[] nonNegativeTipDistance, final double[] nonNegativeForce) {
    ModelFunction modelFunction = new ModelFunction(new MultivariateVectorFunction() {
      @Override
      public double[] value(double[] arg0) {
        Preconditions.checkArgument(arg0.length == 1);
        double contourLength = arg0[0];
        // TODO: Can we reuse this vector?
        double[] value = new double[nonNegativeTipDistance.length];
        for (int j = nonNegativeTipDistance.length; --j >= 0;) {
          value[j] = Wlc.evaluate(kbtOverP, nonNegativeTipDistance[j] / contourLength);
        }
        return value;
      }
    });
    ModelFunctionJacobian jacobian = new ModelFunctionJacobian(new MultivariateMatrixFunction() {
      @Override
      public double[][] value(double[] arg0) {
        Preconditions.checkArgument(arg0.length == 1);
        double contourLength = arg0[0];
        // TODO: Can we reuse this matrix?
        double[][] jacobian = new double[nonNegativeTipDistance.length][1];
        for (int j = nonNegativeTipDistance.length; --j >= 0;) {
          jacobian[j][0] = Wlc.evaluateContourLengthPartialDerivative(
            kbtOverP, nonNegativeTipDistance[j], contourLength);
        }
        return jacobian;
      }
    });
    InitialGuess initialGuess = new InitialGuess(
      new double[] {CommonThings.max(nonNegativeTipDistance) * 1.1});
    Target target = new Target(nonNegativeForce);
    Weight weight = new Weight(CommonThings.nCopies(nonNegativeForce.length, 1.0));
    MaxEval maxEval = new MaxEval(1000);

    final PointVectorValuePair optimum;
    try {
      optimum = estimator.optimize(modelFunction, jacobian, target, weight, initialGuess, maxEval);
    } catch (TooManyEvaluationsException e) {
      logger.log(Level.WARNING, "Could not estimate contour length.", e);
      throw new MyRuntimeException(e.getMessage(), e);
    }
    return optimum.getPointRef()[0];
  }
}
